arXiv:1504.00620v2 [nlin.CD] 21 Apr 2015 


Comparison of Very Smooth Cell-Model Trajectories 
Using Five Symplectic and Two Runge-Kutta Integrators 

Wm. G. Hoover and Carol G. Hoover 
Ruby Valley Research Institute 
Highway Contract 60, Box 601, Ruby Valley 89833, NV USA 
(Dated: April 22, 2015) 

Abstract 

Time-reversible symplectic methods, which are precisely compatible with Liouville’s phase- 
volume-conservation theorem, are often recommended for computational simulations of Hamil¬ 
tonian mechanics. Lack of energy drift is an apparent advantage of such methods. But all 
numerical methods are susceptible to Lyapunov instability, which severely limits the maximum 
time for which chaotic solutions can be “accurate”. The “advantages” of higher-order methods 
are lost rapidly for typical chaotic Hamiltonians. We illustrate these difficulties for a useful 
reproducible test case, the two-dimensional one-particle cell model with specially smooth forces. 
This Hamiltonian problem is chaotic and occurs on a three-dimensional constant-energy shell, 
the minimum dimension for chaos. We benchmark the problem with quadruple-precision trajec¬ 
tories using the fourth-order Candy-Rozmus, fifth-order Runge-Kutta, and eighth-order Schlier- 
Seiter-Teloy integrators. We compare the last, most-accurate particle trajectories to those from 
six double-precision algorithms, four symplectic and two Runge-Kutta. 
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I. INTRODUCTION 


The ongoing computational revolution in physics relies on accurate solutions of funda¬ 
mental equations, Newton’s ( or Lagrange’s or Hamilton’s ) Laws of Motion, in the case of 
classical mechanics. The determinism of these ordinary differential equations is illusory in 
many cases, as typically the equations are “Lyapunov unstable”. Such instabilities grow 
exponentially fast, ~ , where A is the largest Lyapunov exponent of the solution. 

Particle mechanics, our own research interest, provides many examples ranging from 
one-particle chaos to biomolecule simulations using models with many thousands of atomic 
degrees of freedom^. We consider here the simplest particle model for chaos, a one- 
body “cell model” with the periodic four-body cell boundaries shown in Figure 1. The 
resulting motion, approximated with the simplest possible “leapfrog” integrator, described 
below, is generally Lyapunov unstable.-i^ We simplify the initial conditions by starting 
the particle trajectory in the field-free cell interior. We benchmark this problem with 
three quadruple-precision integrators using timesteps chosen to maximize accuracy. We 
compare the resulting benchmark trajectory to six other trajectories from self-starting 
double-precision algorithms typical of molecular dynamics simulations. Five of these 
algorithms are “symplectic”, including the justifiably-popular Leapfrog Algorithm. The 
two others are Runge-Kutta algorithms. 

In the following Sections we describe the specially-smooth differential equations gov¬ 
erning the motion of the wandering cell-model particle, and then quantify the algorithmic 
accuracy with which Leapfrog and the six more sophisticated integrators “solve” this same 
problem. Our conclusions make up the hnal Summary section. 


II. THE CELL MODEL TRAJECTORY IN TWO SPACE DIMENSIONS 

Cell models played a role in models of the liquid state long before the development of 
molecular dynamics.- The geometry treated here is shown at the left in Figure 1. A mass 
point, the “wanderer” particle, moves in a periodic square cell with a motionless fixed 
particle at each of the four vertices. Using periodic boundary conditions the equations of 
motion are : 

X = {px/rn) ; y = {py/m) ] Px = Fx ] Py = Fy . 
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Figure 1; The periodic 2x2 unit cell is shown at the left. The black regions, with potential 
energy greater than one half, are inacessible to the wanderer particle. Initially the wanderer 
is at the origin with velocity (0.6,0.8). Outside the central diamond-shaped region the fixed 
scatterers at the cell corners exert repulsive forces on the wanderer particle. A visually-accurate 
trajectory, calculated with a quadruple-precision fifth-order Runge-Kutta integrator, using five 
million timesteps and dt = 0.00001, is shown at the right with filled circles marking the con¬ 
figurations at times 10, 20, 30, and 40. The open circle corresponds to the maximum time 
t = 50. 



The force on the wanderer is the gradient of the potential function $ , a sum over the 
contributions of the four corner scatterers located at { } : 

4 

<h = 1 — (r — for \r — Til < 1 . 

1 

After advancing the coordinates one timestep dt it is convenient to localize the motion to 
the cell centered on the origin. Whenever the wanderer moves “out”, we replace it “in” 
the basic 2x2 unit cell as follows : 

X < —1 ^ X = X + 2 ] X > +1 ^ X = X — 2 

y<-l^y = y + 2- y>+l^y = y-2. 

We choose initial conditions { x,y,Px,Py } = { 0.0, 0.0, 0.6, 0.8 } and show an accurate 
benchmark solution of the motion equations for a time of 50 at the righthandside of 
Figure 1 . At times of 10, 20, 30, 40, and 50 the benchmark values of {x,y) are : 
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10; +0.321356333887505, +0.585921713605895 
20: +0.81481797353866, -0.572042192203162 
30: -0.040449409487, +0.38290501902 
40: +0.3742439, -0.842854 
50: +0.4696, -0.3568 

Figure 1 shows a unit cell of a periodic two-dimensional lattice in which a single 
particle moves in the held of scattering particles arranged in a hxed square lattice with 
nearest-neighbor spacing of 2. The potential energy maximum of unity is twice the energy 
of the initial condition, shown at the center of the cell. The benchmark solution of the 
motion equations { q = p p = F{q) } is shown at the right. This same accurate 
trajectory was obtained with both the Candy-Rozmus fourth-order and a Runge-Kutta 
hfth-order integrator using 50 million and 500 million timesteps, respectively. The two 
trajectories agree throughout within visual accuracy. At a time of 50 {x,y,px,Py) are : 

{x,y,Px,Py) = (+0.46961,-0.35683,+0.11945,+0.98408) [ CR4 ] ; 

{x,y,Px,Py) = (+0.46962,-0.35682,+0.11948,+0.98408) [ RK5 ] . 

III. SEVEN TYPICAL INTEGRATORS AND THEIR TRAJECTORIES 

We consider seven solution algorithms for the wanderer particle trajectory, [1] Leapfrog 
(symplectic), [2] Fourth-Order Candy-Rozmus Symplectic, [3] Monte Carlo Symplectic, 
[4] Sixth-Order Symplectic, [5] Fourth-Order Runge-Kutta, [6] Fifth-Order Runge-Kutta, 
and [7] Eighth-Order Schlier-Seiter-Teloy Symplectic. For the hrst six of these we use a 
hxed timestep typical of “accurate” molecular dynamics simulations dt = 0.001 . Solutions 
for those six integrators appear in Figures 2-7 . The particle mass is unity and the energy 
<h + K is one half. For the last integrator, which has a trajectory visually identical to that 
of Figure 1 we have chosen timesteps as small as 0.00000001 in order to obtain ten-digit 
accuracy in the wanderer trajectory up to a time of 50 . Let us consider the details of all 
the integrators next. 
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A. Second-Order Time-Reversible Leapfrog Algorithm 


“Symplectic” integrators^^- automatically obey Liouville’s Theorem by advancing the 
solution of Hamiltonian problems in time according to a series of phase-volume-conserving 
shears. Symplectic algorithms alternate steps advancing the coordinates and momenta in 
time. The simplest example is equivalent to the Stormer-Verlet “leapfrog algorithm” 

{ q = q + (p * dt/2) ; p = p + (F/m) * dt ; q = q + (p * dt/2) } <(—)■ 

{ qn+i - 2qn + qn-i = {F/m)n } . 

This algorithm is said to be “second order”—, with a hxed-time coordinate error of order 
tdF for t « 271/dF when applied to the simple harmonic oscillator. It is time reversible 
in that changing +dt —)■ —dt gives the same trajectory points either forward or backward 
in time. 

How does the simulation begin? Starting out at the origin, with the wanderer speed 
equal to unity and a hxed timestep dt = 0.001, the hrst 420 steps leave the momenta 
unchanged and becomes 1.0004. During the 421st step the upper right scatterer is 
contacted and begins to repel the wandering particle with a force : 

= 8(x - 1)(1 - ; Fy = 8(y - 1)(1 - ; r = (1 - x,l -y) , 

where x and y are the wanderer coordinates. 

After an elapsed time t we reverse the sign of the time so as to integrate backward to 
see how closely the wanderer returns to its initial location. So long as f < 47 we hnd 
that the trajectory reverses to within a distance 0.01 of the origin. We will see that this 
retracing of steps does not guarantee a match with the accurate trajectory shown at the 
right in Figure 1. Both the trajectory reversal and the conservation of energy are poor 
diagnostics for trajectory accuracy, where accuracy means reproducing correct values of 
the coordinates x(t),y(t) . 

B. Fourth-Order Time-Reversible Symplectic Integrator 

Higher-order algorithms, with hxed-time integration errors of order dt^, dt^, dt^ ... can 
be developed from Taylor’s series about t giving small increments in the coordinates and 
momenta as three-, four-, hve- ... term series in dt . Candy and Rozmus’ fourth-order 
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Figure 2: The Leapfrog integrator reproduces the accurate x and y coordinates within 0.01 
for an integration time of 18. The energy at that point ( where the trajectory color changes, 
indicated by a star ) is in error in the seventh decimal place. Here and elsewhere the cited double¬ 
precision times are truncated to integers because different implementations, such as varying the 
order of the operations, could change these numbers. 

integrator (with an error of order dt^ at a fixed not-too-large, time) is a simple example, 
cited in the very useful summary paper by Gray, Noid, and Sumpter- : 

q = q -|- 0.6756036p * dt ; p = p + 1.3512072(F/m) * dt ; 

q = q — 0.1756036p * dt ; p = p — 1.7024144(F/m) * dt ; q = q — 0.1756036p * dt ; 

p = p -F 1.3512072(F/ni) * dt ; q = q + 0.6756036p * dt . 

Reference 7 gives the analytic forms of all of the coefficients. Notice that the coefficients 
incrementing the coordinates sum to unity as do also those incrementing the momenta. 
Each timestep requires three separate evaluations of the forces. 

C. Monte-Carlo Time-Reversible Symplectic Integrator 

Although it is usual to provide coefficients in integration algorithms to many signih- 
cant hgures, in most cases an approximate rendition is sufficient. It is quite possible to 
develop algorithms with a Monte Carlo method, adjusting the coefficients to minimize 
the trajectory error for the simple harmonic oscillator problem. An integrator requiring 
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Figure 3: The Candy-Rozmus fourth-order symplectic trajectory exhibits a color change at a 
time of 34, the maximum for which the coordinate errors are less than 0.01 . The energy error 
at that time is in the twelfth decimal place. The maximum time at which a reversed trajectory 
returns to the origin within 0.01 is t = 42 . 

five force evaluations per timestep was developed by Monte Carlo sampling^ adjusting the 
coefficients subject to the constraints of time reversibility and normalization so that the 
Monte Carlo trajectory optimization occurs in a four-dimensional space. The resulting 
integrator was successful in modelling many-body dynamics but is here applied to the 
cell-model problem of Figure 1 : 

q = q -|- 0.005904p * dt ; p = p + 0.171669(F/m) * dt ; 

q = q -f 0.515669p * dt ; p = p — 0.516595(F/ni) =t= dt ; 
q = q — 0.021573p * dt ; p = p + 1.689852(F/m) * dt ; q = q — 0.021573p * dt ; 
p = p — 0.516595(F/ni) =t= dt ; q = q + 0.515669p * dt ; 
p = p -(- 0.171669(F/m) * dt ; q = q + 0.005904p * dt . 

The cell model trajectory using this Monte Carlo integrator is illustrated in Figure 4. 

D. Yoshida’s Sixth-Order Time-Reversible Integrator 

Yoshida developed and applied a general technique for hnding a variety of higher- 
order symplectic integrators.— His sixth-order time-reversible integrator advances the 
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Figure 4: The color change in the trajectory from the Monte Carlo symplectic integrator occurs 
at a time of 31, after which the coordinate errors exceed 0.01. The energy error there is in the 
thirteenth digit. For this integrator a trajectory reversed at a time of 43 will return to the origin 
with coordinates recurring within 0.01. 

coordinates (Ag oc pdt) eight times per timestep, using the symmetric (so as to guarantee 
time-reversibility) set of eight coefficients which sum to unity : 

+0.39225680523878, +0.51004341191846, -0.47105338540976, +0.06875316825252, 

+0.06875316825252, -0.47105338540976, +0.51004341191846, +0.39225680523878. 

Between the successive coordinate updates there is a force calculation and an update of 
the momenta (Ap oc Fdt), using seven coefficients, which likewise sum to unity : 

+0.78451361047756, +0.23557321335936, -1.1776799841789, +1.3151863206839, 

-1.1776799841789, +0.23557321335936, +0.78451361047756. 


E. Fourth-Order and Fifth-Order Runge-Kutta Integrators 

Runge-Kutta integrators ( circa 1900, as described in Wikipedia ) advance both co¬ 
ordinates and momenta simultaneously in a series of stages within each timestep dt. As 
in the symplectic case the variables at time t + dt are expressed as series in dt, putting 
conditions on the summed-up coefficients for each power of dt to be treated correctly by 
the algorithm. 






Figure 5: This double-precision trajectory is based on Yoshida’s time-reversible sixth-order 
integrator with a timestep dt = 0.001 . There is a color change at t = 36 , indicating the 
degradation of trajectory accuracy to ±0.01 despite the negligible energy error in the fourteenth 
decimal place. Trajectory reversal at a time of 42 returns to the origin within coordinate errors 
of 0.01. 

The main advantage of Runge-Kutta methods is that they can be applied to arbi¬ 
trary sets of ordinary differential equations, not just those from Hamiltonian mechanics. 
The fourth-order “classic” Runge-Kutta method has been a standard workhorse model 
for solving sets of coupled ordinary differential equations for 100 years. Applied to the 
harmonic oscillator the fourth-order algorithm suffers a loss in energy proportional to 
the fifth power of the timestep. The hfth-order Runge-Kutta integrator behaves in the 
opposite manner with the energy increasing rather than decreasing. 

Hybrid “adaptive” models, incorporating both fourth- and fifth-order algorithms, pro¬ 
vide a simple means for the automatic control of integration errors. The harmonic oscilla¬ 
tor is an excellent test case of integrator accuracy where Lyapunov instability is absent.— 
Figure 6 illustrates the same cell-model orbit for the classic fourth-order Runge-Kutta 
integrator. Figure 7 shows a fifth-order Runge-Kutta integrator : 
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Figure 6: The fourth-order Runge-Kutta trajectory using double precision and a timestep 
dt = 0.001 provides coordinates accurate within 0.01 through a time of 35, indicated by the 
color change at the star. The energy error at that point, 10“^^, is negligible. Changing the sign 
of the timestep at t = 42 , +dt —dt , returns the trajectory to the origin within a precision 
of 0.01 . 
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Figure 7: The fifth-order Runge-Kutta trajectory using double precision and a timestep dt = 
0.001 provides coordinates accurate within 0.01 through a time of 37, indicated by the color 
change at the star. The energy error at that point is 10“^^^ . For times less than 42 reversing the 
trajectory, by setting +dt —>■ —dt , returns the trajectory to the origin with precision 0.01. This 
integrator is the best of the double-precision integrators tested here. Any one of the five higher- 
order integrators is accurate for about twice the time of the second-order Leapfrog integrator. 
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ypl = yp[ y ; 


yp2 = yp[ y+ (dt/2) * ypl ] 

yp3 = yp[ y + (dt/16) * (3ypl + yp2) ] 

yp4 = yp[ y+ (dt/2) * yp3 ] 

yp5 = yp[ y + (dt/l6) * (-3yp2 + 6yp3 + 9yp4) ] 

yp6 = yp[ y + (dt/7) * (ypl + 4yp2 + 6yp3 - 12yp4 + 8yp5) ] 

y = y + (dt/90) * (7ypl + 32yp3 + 12yp4 + 32yp5 + 7yp6) 

Here yp [ ... ] represents the righthandside of the vector differential equation y = y' 

where the six force evaluations in each timestep are indicated by { ypl, yp2, ... yp6 } . 

Both Runge-Kutta integrators return to the origin with errors no more than 0.01 with 
reversal at time 42. Forward in time their trajectories are accurate through times of 35 
and 37, the last being the best of the double-precision integrators. The energy errors for 
the two Runge-Kutta integerators are in the thirteenth and fourteenth decimal places. 

F. An Eighth-Order Time-Reversible Symplectic Integrator 

Ernst Teloy, Christoph Schlier, and Ansgar Seiter developed and implemented a use¬ 
ful eighth-order time-reversible symplectic integrator with 17 force evaluations per step. 
Applied to the harmonic oscillator the rms coordinate error increases by about eight or¬ 
ders of magnitude when the timestep is increased by a factor of ten, consistent with an 
eighth-order method. 

For the reader’s convenience we reproduce here the 18 coefficients required to imple¬ 
ment the method. They can be found quoted to 35 decimal places at Christoph Schlier’s 
Freiburg website or in Reference 12 . This precision is steadily reduced, digit by digit, 
through Lyapunov instability, described in more detail in Section IV. In the cell-model 
case the rate of precision loss is 0.7 , one binary bit per unit time. Accordingly, for the 
eighth-order integrator in quadruple precision at time 50 we would expect an exponen¬ 
tially amplihed error of order 10’ 32 ^ 2^0 ~ 10 . In fact, we hud a trajectory error of 

order 10“^° using a timestep of 10“® , as is shown below. 

Even so the eighth-order integrator with dt = 0.001 loses only seven of the original 
35 digits in energy along with twenty digits in position when run forward and backward 
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for 50,000 steps to match the time illustrated in all the Figures. As was illustrated and 
emphasized in References 12 and 13 energy conservation and trajectory reversibility are 
both of them misleading diagnostics of trajectory accnracy. It is only throngh a study of 
convergence that trajectories can be validated. For the eighth-order symplectic integrator 
the timestep dependence of the {x, y) coordinates at time 50 is as follows : 


dt = 0.00100000 ^ (0.48704 51729, +0.13435 10401) 


dt 

dt 

dt 

dt 

dt 


0.00010000 ^ (0.48185 80396, 
0.00001000 ^ (0.46961 32018, 
0.00000100 ^ (0.46961 40145, 
0.00000010 ^ (0.46961 40143, 
0.00000001 ^ (0.46961 40142, 


-0.32559 07485) 
-0.35683 11339) 
-0.35682 95856) 
-0.35682 95861) 
-0.35682 95862) 


For the convenience of the reader we reproduce the integrator coefficients here from Ref¬ 
erence 12, together with a short harmonic-oscillator program to demonstrate their use. 
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G. Schlier-Seiter-Teloy Integrator Coefficients 


c( 1) 

= 

+0.04463 

79505 

23590 

22755 

91399 

96257 

33590 

dOO 

c( 2) 

= 

+0.13593 

25807 

16909 

59145 

54326 

42134 

95574 

dOO 

c( 3) 

= 

+0.21988 

44042 

71470 

72254 

44553 

50696 

06167 

dOO 

c( 4) 

= 

+0.13024 

94678 

05238 

28601 

62119 

37781 

96846 

dOO 

c( 5) 

= 

+0.10250 

36569 

39750 

69608 

26124 

10077 

79814 

dOO 

c( 6) 

= 

+0.43234 

52186 

93585 

47487 

98325 

78848 

77035 

dOO 

c( 7) 

= 

-0.00477 

48291 

69168 

81658 

02248 

90639 

62934 

dOO 

c( 8) 

= 

-0.58253 

47690 

40408 

45493 

11283 

79308 

61212 

dOO 

c( 9) 

= 

-0.03886 

26428 

21118 

17697 

73742 

08751 

89743 

dOO 

c(10) 

= 

+0.31548 

72853 

79404 

79698 

27360 

37972 

74199 

dOO 

c(ll) 

= 

+0.18681 

58374 

32971 

55471 

52615 

35039 

72746 

dOO 

c(12) 

= 

+0.26500 

27549 

90620 

83398 

34600 

29630 

79872 

dOO 

c(13) 

= 

-0.02405 

08473 

57473 

61993 

57358 

79824 

07554 

dOO 

c(14) 

= 

-0.45040 

49249 

97722 

51180 

92289 

67121 

51891 

dOO 

c(15) 

= 

-0.05897 

43301 

55923 

86914 

57532 

39267 

66330 

dOO 

c(16) 

= 

-0.02168 

47617 

18613 

35324 

93438 

86847 

07580 

dOO 

c(17) 

= 

+0.07282 

08003 

35901 

28173 

76189 

26412 

34244 

dOO 

c(18) 

= 

+0.55121 

42963 

41970 

67334 

40560 

13815 

94315 

dOO 


Oscillator program with q,p,dt and 18 c(i) ; 
do i = 1,17,2 
q = q + c(i)*p*dt 
p = p - c(i+l)*q*dt 
enddo 

do i = 17,1,-2 
q = q + c(i)*p*dt 
if(i.gt.l) p = p - c(i-l)*q*dt 
enddo 
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IV. LYAPUNOV INSTABILITY IN THE CELL MODEL 


For chaotic systems the algorithmic accuracy of numerical integrators deteriorates ex¬ 
ponentially rather than linearly in the time^. The underlying exponential Lyapunov in¬ 
stability of dynamical systems is easily measured by following the motion of a “reference” 
trajectory in the usual way, for instance with any one of the seven algorithms discussed 
here. An additional “satellite” trajectory, separated from the reference by a small length 
5o , is also followed using the same algorithm. At the end of each timestep the separation 
is rescaled, maintaining the length of the offset between the trajectories constant, but 
allowing the direction to vary : 

6{t + dt) = [ Tsit -|- dt) — Trit + dt) ] ; 

rg —>■ rr{t + dt) + 6{t + dt)[ 5 q/\ 6{t + dt) \ ] . 

The largest Lyapunov exponent is simply the average value of the growth rates measured 
at the ends of every timestep prior to rescaling : 

Ai = ( (l/dt) ln[ I 6{t + dt) |/^o ] ) • 

Previous studies of this cell model,- with the same initial condition, have shown that the 
largest Lyapunov exponent is about 0.7. This means that an error of the order 10“^® at 
the initiation of a run of length 50 will increase by a factor of ~ 10^^ . 

This exponential growth rate explains why it is that all of the double-precision inte¬ 
grators fail, from the standpoint of reproducing a reversible trajectory, at about the same 
time, at about half the time where quadruple-precision trajectories fail. It is because 
these trajectories are just approximations that the most sophisticated biomolecule simu¬ 
lations are based on the rudimentary leapfrog algorithm rather than more sophisticated 
algorithms. 

Of course, even the slightest difference in the error prior to amplihcation will yield 
a different history. Just summing the particle interactions in a different order leads to 
qualitatively different histories once the Lyapunov instability rises to the level of visibility, 
an increase of 16 digits for routine double-precision simulations. The phase-shift errors 
in all of the algorithms discussed here can be measured by choosing the initial velocity 
[^1/2, ^Jl/2) for which the roundoff errors in the x and y directions are identical. 
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If high accuracy is required, as in astronomical simulations, multiple precision can be 
employed, as demonstrated by Lorenz Attractor simulations using a precision of thousands 
of decimal digits. But, as Joseph Ford was fond of pointing out, Lyapunov instability 
is incompatible with high accuracy. Doubling the number of signihcant hgures in the 
integration algorithm only doubles the time for which the simulation is accurate. 

Recently Hanno Rein and David Siegel^ developed and implemented a relatively com¬ 
plicated fifteenth-order integrator for gravitational problems with the provocative title 
“A Fast, Adaptive, High-Order Integrator for Gravitational Dynamics, Accurate to Ma¬ 
chine Precision Over a Billion Orbits”. Evidently this integrator is not at all intended for 
long-time applications to chaotic problems, where errors grow exponentially with time. 
Conversations with Ben Leimkuhler and Mark Tuckerman, both of whom summarily dis¬ 
miss the use of Runge-Kutta techniques, due to their monotonic energy drift, plus the 
appearance of Rein and Siegel’s high-order long-time work led to the present article. 

V. RECENT DEVELOPMENTS AND SUMMARY 

To summarize, for simple chaotic simulations ( such as classical fluids ) symplectic inte¬ 
grators attain accuracies similar to those obtained with Runge-Kutta integration and are 
primarily limited by Lyapunov instability. Although energy conservation and trajectory 
reversibility characterize symplectic integrators, those properties do not ensure trajectory 
accuracy. The reversibility of the double-precision leapfrog integrator, to a time of 47 and 
back, exceeds that of all the more accurate double-precision integrators. 

For us it was illuminating to find that the humble Leapfrog integrator, presumably 
nearing its 330th anniversary^, is nearly as useful as are its more complex relatives, and 
is certainly far more economical. For higher accuracy there is little distinction between 
the symplectic and the Runge-Kutta integrators for chaotic problems, because both types 
lose accuracy at the very same rate, determined by the maximum Lyapunov exponent. 

It is significant that all of the integrators used here conserve energy almost perfectly 
for the benchmark problem. They also reverse back to the initial conditions even when 
their trajectories are inaccurate. One takeaway message from these simulations is the one 
to which Joseph Ford devoted much thought and many thought-provoking words, among 
them these taken from Reference 16 : 

“Newtonian determinism assures us that chaotic orbits exist and are unique. 
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but they are nevertheless so complex that they are humanly indistinguishable 
from realisations of truly random processes.” 

Liao has confronted the Lyapunov instability problem headon for the Lorenz 
Attractor.— By using 3500-term series expansions coupled with 4180-digit arithmetic he 
followed the evolution of the Lorenz Model to a time of 10,000. Like the continuing 
discovery of the digits of vr this activity will last as long as mankind. 

Lyapunov instability often shows up in peculiar places. Simply changing the order of 
operations in adding up forces or in computing the weights of contributions to differential 
equations’ righthandsides can provide the seeds from which macroscopic change develops. 
We learned this lesson in simulating the collisions of mirror-image manybody drops and 
crystals. To retain accurate mirror symmetry it was necessary to symmetrize the force 
calculations at every timestep.— 
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